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SUMMARY 


The results of an investigation of boundary layers close to the stagnation point 
of an oscillating airfoil are reported. Two procedures for generating initial condi- 
tions — the characteristic-box scheme and a quasi-static approach — were investigated, 
and the quasi-static approach was shown to be appropriate provided the initial region 
was far from any flow separation. With initial conditions generated in this way, the 
unsteady boundary- layer equations were solved for the flow in the leading-edge region 
of a NACA 0012 airfoil oscillating from 0° to 5°. Results were obtained for both 
laminar and turbulent flow, and, in the latter case, the effect of transition was 
assessed by specifying its occurrence at different locations. The results demonstrate 
the validity of the numerical scheme and suggest that the procedures should be applied 
to the calculation of the entire flow around oscillating airfoils. 


INTRODUCTION 


The calculation of boundary-layer characteristics of an oscillating airfoil dif- 
fers from the usual unsteady flow calculations in that difficulties are caused by the 
translation of the stagnation point in space and time. In particular, it is essential 
to develop a procedure to generate initial conditions in the immediate vicinity of 
the moving stagnation point and to account for the flow reversal that occurs in this 
region. 

The study reported here is the continuation of the work described in reference 1. 
It is one phase of a study that will be extended later to compute the complete bound- 
ary layer and inviscid flow characteristics of an oscillating airfoil in order to 
improve understanding of the dynamic-stall problem. In the present study, we focus 
our attention on the calculation of boundary layers near the stagnation point of an 
oscillating airfoil. We consider both laminar and turbulent flows and two different 
procedures for generating the initial conditions in the (t,y) plane. 

The following section describes the basic equations, turbulence model, the 
initial conditions and the solution method. The details of the numerical procedure 
were discussed in reference 1 and are not repeated here. 
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The third section presents the results for NACA 0012 airfoil. Calculations were 
first performed with two procedures to investigate the prediction of initial condi- 
tions in the (t,y) plane. They were limited to laminar flow and to the neighborhood 
of the leading edge of the airfoil. The next set of calculations involved the 
boundary-layer behavior of the NACA 0012 airfoil oscillating between angles of attack 
of 0° and 5° for laminar, transitional, and turbulent flows and for a chord Reynolds 
number of 3><10 6 . The transition location was varied in order to investigate its 
influence on flow separation. 


BASIC EQUATIONS 


Boundary-Layer Equations 


The boundary-layer equations for an incompressible laminar or turbulent flow on 
an oscillating airfoil are well known and, with the eddy-viscosity (e m ) concept, can 
be written as 
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where x denotes distance along the surface of the airfoil, y is distance along the 
normal, and b = v + e m . In the absence of mass transfer, equations (1) and (2) are 
subject to boundary conditions given by 


y = 0 ; u = v = 0 ; y = 6 ; u= u e (x,t) (3) 

The presence of the eddy viscosity e m requires a turbulence model; we use the alge- 
braic eddy-viscosity formulation developed by Cebeci and Smith (ref. 2). According 
to this formulation, e m is defined by two separate formulas. In the inner region of 
the boundary layer (e m )i is defined as 
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In equation (4), y tr is an intermittency factor that accounts for the transi- 
tional region that exists between a laminar and turbulent flow. It is defined by 
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Here x tr is the location of the start of transition and the empirical factor G, 
which has the dimensions of velocity/ (length) 2 , is given by (ref. 2) 
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The transition Reynolds number is defined as ^x tr = (u e x/v) tr . 
In the outer region, (e^o defined as 


(e ) = 0.0168 

mo 



( 8 ) 


The boundary between the inner and outer regions, y c , is established by the continuity 
of the eddy-viscosity formulas. 


If initial conditions in the (t,y) plane are given at a station x Q on the 
upper surface of the airfoil and satisfy the condition u > 0 and, in addition, 
initial conditions are given in the (x,y) plane at t = 0 , then the solutions of 
equations ( 1 ), (2), and (3) may be integrated in x > x Q until they break down (flow 
separation). A similar remark applies to the lower surface except that u < 0. The 
initial conditions at t = 0 can be generated for both surfaces if steady conditions 
are assumed to prevail at that time. It is only necessary to solve the appropriate 
equations which, in this case, are given by equation ( 1 ) and by 


There is no problem with the initial conditions for equations (1) and (9) since the 
calculations start at the stagnation point x = x g , where u e and u are zero for 
all y . 

Unlike steady flows, where u e and u are zero for all y at the stagnation 
point, the stagnation point is not fixed in an unsteady flow; although u e is zero, 
we cannot assume a priori that u is also zero. We may avoid these difficulties by 
using an implicit method, but now we are faced with the problem of generating a 
starting profile on the new time-line. 

A convenient and accurate procedure for calculating the first velocity profile 
at the new time-line has been developed (ref. 1 ); it involves the use of the 
"characteristic-box 11 scheme developed by Cebeci and Stewartson (in ref. 3) and 
described in ref. 1. Another procedure, though not as accurate, is to use a quasi- 
steady approach in the immediate vicinity of the stagnation point region. We shall 
discuss both procedures later (Results and Discussion section). 


As in previous studies (see, for example, ref. 4), we use similarity variables 
to transform the governing equations before we seek their solution. For a steady 
flow, we use the Falkner-Skan transformation defined by 


Initial Conditions 
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Transformed Equations 
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where \p is the usual definition of stream function that satisfies the continuity 
equation (1), that is 
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With this transformation, equation (9) and its boundary conditions, equation (3), can 
be written as 
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where primes denote differentiation with respect to n , and m denotes a dimension- 
less pressure-gradient parameter defined by 


m = 
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For unsteady flows, we use a transformation similar to that defined by equa- 
tion (10) except that u e is now a function of both x and t, and the dimensionless 
stream function F is a function of x, t, and we let 


[u (x,t)] l/2 

C = — y > V = [u (x,t)vx] 1 / 2 F(x,t,c) 
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With this transformation, it can be shown that the continuity and momentum equations 
and their boundary conditions for unsteady incompressible flows can be written as 
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Here primes now denote differentiation with respect to £ and 
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Solution Procedure 

We use Keller’s box method to solve the governing equations of the previous 
section. This is a two-point finite-difference method which has been used to solve a 
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wide range of parabolic partial-differential equations, as discussed in reference 4. 
The solution procedure for equations (12) and (13) is identical to that described in 
reference 5. The solution procedure employing the characteristic box scheme to 
generate the first velocity profile at a new time-line is described in reference 1. 
For unsteady flows, where we now solve equations (16) and (17), we use the solution 
procedure described in reference 6. In regions where there are no flow reversals 
across the layer, we use the "standard box*' scheme and in regions where there is flow 
reversal, we use the "zig-zag" scheme. 


A Model for External Velocity Distribution 


The solution of boundary-layer equations requires that the external velocity 
distribution be specified. Since the present effort is directed toward solutions 
near the leading edge of the airfoil, a local model for the potential flow has been 
chosen in the place of a full-potential flow code. We first consider an ellipse with 
major axis 2a and thickness ratio T, where 2aX and 2AY, respectively, measure 
distance along and perpendicular to the major axis from one apse, that is, from the 
nose. The equation of the ellipse is then 
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and the velocity distribution near the nose is given by (ref. 6) 

= + (1 + T) (X 1 ^ 2 t ot) 

" [X + ( 1 /4)T 2 ] 1/2 


( 20 ) 


where a is proportional to the angle of attack in radians. Here T denotes the 
thickness ratio (E b/a) and (+) denotes the upper surface and (-) the lower surface. 
We note that equation (20) is valid only when a = 0(T), X = 0(T 2 ), and the location 
of the stagnation point X s is X s = a 2 . 

The external velocity distribution of a symmetrical airfoil in the neighborhood 
of the nose can also be represented by an expression similar to equation (20). It 
is only necessary to let T denote the thickness ratio of the equivalent ellipse and 
replace T 2 /4 by R/2, where R is the nose radius. For example, the nose region 
of a NACA 0012 airfoil at an angle of attack of 5° can be represented by the ellipse 
by matching three points (0, 0.2, and 0.12); then the appropriate form of equa- 
tion (20) is 
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Figure 1 shows a comparison of the external velocity distribution computed by the 
inviscid flow theory and that computed by equation (21) for the lower and upper sur- 
faces of a NACA 0012 airfoil. Equation (21) is thus a satisfactory fit for the 
velocity distribution extending to approximately 20% chord. 

Equation (21) can also be used to approximate the external velocity distribution, 
when the angle of attack is varying sinusoidally, by writing it as 
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( 22 ) 


— = +1 2065 [ xl/2 * 0-957a(l + A sin a>t)] 

U °° (X + 0.00613) 1 / 2 

Here physical time is u^t/c, the period of oscillation of the airfoil is 
2 7T c / U qqW , c (= 2a) is the chord, and A is a free parameter. We note from equa- 
tion (8) that at time t = 0, the external velocity distribution corresponds to that 
given by steady state. By choosing A = 1 and a = 5°, we can compute the velocity 
distribution at angles of attack ranging from 0° to 10°. The term 

0.957a(l + A sin ait) (23) 

can be interpreted as an effective angle of attack, a e ff(t). With this definition, 
equation (23) becomes 
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Figure 2 shows a comparison of the external velocity distribution computed by Neumann 
and by equation (24) for a = 0° and 10°. As can be seen, equation (24) is a satis- 
factory representation of the velocity distribution, especially near the leading 
edge. We note from the results that the agreement begins to deteriorate on the upper 
surface as a e ff begins to increase. However, equation (24) is a convenient formula 
for generating the external velocity distribution in the immediate vicinity of the 
stagnation point and for testing the computer program. 


RESULTS AND DISCUSSION 


Calculation of Initial Conditions by Two Separate Procedures 


The procedure, based on the so-called characteristic-box scheme and developed 
for calculating the initial conditions in the (t,y) plane, was tested for a model 
problem (ref. 1) in which the external velocity distribution was given by 


£ + 1 + A sin cot) 

(1 + 5) a ] 1/2 


(25) 


Here X = (1/4)T 2 £ 2 , £ = 0(1), and A and denote parameters that need to be speci- 

fied. The latter can be regarded as a reduced angle of attack; if it is increased 
beyond 1.155, the steady-state solution at t = 0 separates on the upper surface of 
the airfoil. Calculations were carried out with this procedure (ref. 1) for a test 
case — K 0 = 0.10, A = 1, co = tt/ 4 — for a limited range of x(|x| <0.3), and the formal 
validity and efficacy of the numerical scheme was established. 


When this procedure was used for the external velocity distribution given by 
equation (24) with co = 0.1, a value typical of conditions of dynamic stall, the com- 
puting time of the numerical procedure increased significantly and the unsteady 
effects on the stagnation-point behavior decreased. For this reason, an alternative 
procedure based on a quasi-steady approach was developed and its results were com- 
pared with the predictions of the procedure of reference 1. With the quasi-steady 
approach, the stagnation point was computed from equation (24) and the solution of 
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the boundary-layer equations was obtained up to X = 0.005 (= X Q /c). Because of the 
rapid variation of the external velocity distribution in the stagnation-point region, 
extremely small values in Ax corresponding to 0.001 were taken. Since the quasi- 
steady approach is not valid for flows approaching separation, x Q /c was chosen to 
be just upstream of the pressure minimum and the calculations were then limited to 
the region where the pressure gradient was favorable. 

Figure 3 shows a comparison of the predictions of both procedures. The calcula- 
tions were performed with initial conditions in the (x,y) plane corresponding to a 
steady flow with an effective angle of attack a of 5° and for the external velocity 
distribution given by equation (24). The parameter A was set equal to 1 and the 
reduced frequency was equal to 0.1. The largest effective angle of attack was 10°, 
corresponding to rnt = tt/2, and the calculations were limited to the quarter cycle 
from zero to tt/2. It is expected that any differences between the results of the 
two procedures would be greatest in this range, especially close to tt/2. 

The results shown in figure 3 indicate that the computed values of local skin- 
friction coefficient obtained by both procedures are in remarkable agreement with 
each other. (The small difference in Cf stems from the use of different 
Ax-spacings . ) The procedure that uses the characteristic-box approach requires much 
smaller Ax-spacings than the quasi-steady approach, without offering significant 
improvement in prediction accuracy. The quasi-steady approach is easier to use (than 
the characteristic-box approach), less demanding of computer time, and of equivalent 
accuracy in the region of the stagnation point; as a result, we found it convenient 
for use in the present study. This conclusion is, however, limited to the generation 
of initial conditions near the stagnation point and to situations far from separation. 

It is of interest to investigate the reason for the success of the quasi-steady 
approximation in this study in contrast to the earlier one (ref. 1), in which the 
behavior of the stagnation line is too complicated to be satisfactorily approximated 
in this way. Examination of equation (24) shows that the range of values of X over 
which u e varies by a significant amount is quite small; for a e ff =5°, u e 
increases from zero at X 1 ' 2 = -0.084 to a maximum of l.77u O0 at X = 0.005, the 
corresponding values of x/c being -0.082 and +0.073. Thus, in equation (2), 
u(9u/8x) " 20u^/c, and 3u/8t ~ (u)u 00 /c)u 00 = 0.1u^/c. Hence, in retrospect it is not 
surprising that the unsteady term in equation (1) may be neglected in comparison with 
the steady inviscid terms when computing the boundary-layer structure near the for- 
ward stagnation point. Further downstream, u(3u/8x) becomes much smaller and the 
neglect of 3u/8t is not necessarily justifiable. Indeed, if separation occurs, it 
is crucial that the 3u/3t term be retained in the integration; otherwise, the 
solution will develop a singularity when the skin-friction vanishes. 

The results of the present study may be generalized into the working rule that if 
the frequency of oscillation is u)*/2tt, if R is the nose radius of the airfoil, and 
if co^R/uoo <0.1, then the quasi-steady approximation may be used to compute the 
boundary layer near the nose on the compression side, and as far as the pressure 
minimum on the suction side. For the calculation reported in reference 1, the value 
of a) in equation (25) was taken to be tt/4. Thus the two inertial terms, 3u/8t 
and u(3u/3x), in equation (1) are approximately of the same order of magnitude near 
the forward stagnation point and neither could safely be neglected. 
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Boundary-Layer Behavior of a NACA 0012 Airfoil Near the Leading Edge 

With the initial conditions computed by the quasi-steady approach, calculations 
were performed to investigate the boundary- layer behavior of the NACA 0012 airfoil 
near its leading edge at different angles of attack for a reduced frequency of 0.1 
and for a chord Reynolds number of 3><10 6 . Calculations were first performed for 
angles of attack in the range of 0° to 5°. The expression for the effective angle 
given by equation (23) was written as 

0.957a ,, N 

a eff = 2 (1 ” cos wt ) (26) 

so that the steady-state calculations start for zero angle of attack and, for half a 
cycle, unsteady flow calculations were performed as the angle of attack was varied 
from 0° to 5° by taking a = 5° in equation (26). Results were obtained for both 
laminar and turbulent flows, with the transition point specified at (X/c) tr = 0.06. 

Figure 4 shows the variation of wall-shear parameter and displacement 

thickness 6*, defined by 



with X for three effective angles: a e ff = 0, 2.5, and 5. We see that the location 

of minimum wall shear occurs at the transition point, namely (X/c) tr = 0.06, which 
moves upstream with increasing angle of attack, for example, to X/c = 0.05 at 
a e ff =5°. The displacement thickness reaches a maximum, for example, at X/c = 0.06 
for a e ff = 0°, reduces to a minimum after transition, and increases again with 
further increase in X. The effect of angle of attack on displacement thickness is 
pronounced, and at a e ff = 5°, the displacement thickness reaches a maximum of 
X/c = 0.06, with a large dip in the AX/c range of 0.03, before it increases again. 

Calculations were next performed for the same angle-of-attack range, but this 
time the effective angle given by equation (23) was written as 

a 0ff = ^ + sln wt ) ( 28 ) 

so that the steady-state calculations started at an angle of attack of 2.5°. Unsteady 
calculations were performed for one cycle by taking 10° increments in time. The 
transition point was specified at two different X-locations to determine its effect. 

Figures 5 and 6 show the variation of wall-shear parameter F^ and displacement 
thickness 6* with X for the same effective angles as those in figure 4, but with 
(X/c) tr = 0.04 and 0.07, respectively. The minimum value of wall shear in figure 5 
occurs at the transition point, namely (X/c) tr = 0.04 for steady-state conditions 
with a e ff =2.5° and moves downstream with decreasing angle of attack; for example, 
to X/c = 0.044 at a e ff = 0° and upstream with increasing angle of attack to 
X/c = 0.036 at a e ff =5°. As in figure 4, the effect of angle of attack on dis- 
placement thickness is pronounced as a e ff increases from 0° to 5° causing a large 
dip in 6*. 
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The results of figure 6 were obtained for a delayed transition location of 
(X/c)tr = 0.07 and show similar behavior of wall shear and displacement thickness 
with X at different angles of attack. Since the region of laminar flow has 
increased, the values of wall shear are lower, indicating that the flow may approach 
separation conditions with further delay of transition. The minimum value of wall 
shear again occurs at the transition point, (X/c) tr = 0.07, with increasing or 
decreasing angle of attack as it moves upstream, for example, to X/c = 0.065 at 
a eff = 0° and to X/c = 0.06 at ot e ff = 5°. 

CONCLUDING REMARKS 


A method for calculating the behavior of laminar and turbulent boundary layers 
on an oscillating airfoil has been developed and used to obtain results for the 
boundary layer around a NACA 0012 airfoil oscillating between 0° and 5°. Two proce- 
dures for generating the initial conditions in the (t,y) plane were investigated — 
a characteristic box scheme and a quasi-static approach. The quasi-static approach 
was shown to be preferable provided the X-location was far from flow separation. 

The boundary-layer results were obtained by solving the unsteady-flow equations for 
different angles of attack for both laminar and turbulent flows. They allow compari- 
son of laminar and turbulent flow and for the latter, quantify the effect of changing 
the location of transition. They are presented in terms of displacement thickness 
and wall-shear parameter, both of which show large differences at the larger angles 
of attack. 
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Figure 1.- Comparison of inviscid velocity distribution computed by Neumann (solid 
line) and those computed by equation (21) (dashed lines) for a = 5 (NACA 0012 
airfoil) . 
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Figure 2.- Comparison of inviscid velocity distribution computed by Neumann (solid 
lines) and those computed by equation (24) (dashed lines) for a = 0°, 10°. 
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Figure 3.- Comparison of the predictions of two separate procedures for computing 
initial conditions in the (t,y) plane. Solid lines denote the solutions obtained 
by the characteristic box scheme described in reference 1; the dashed lines 
denote those obtained with the quasi-steady approach. 




(a) Wall-shear parameter F*. (b) Displacement thickness 6*. 

Figure 4.- Variation of wall-shear parameter and displacement thickness with X/c for 
NACA 0012 airfoil with transition set at X/c = 0.06; the effective angle of 
attack is given by equation (26) . 
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(a) Wall-shear parameter F”. (b) Displacement thickness 6*. 

Figure 5.- Variation of wall-shear parameter and displacement thickness with X/c 
for NACA 0012 airfoil with transition set at X/c = 0.04; the effective angle of 
attack is given by equation (28). 
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Figure 6.- Variation of wall-shear parameter and displacement thickness with X/c for 
NACA 0012 airfoil with transition set at X/c = 0.07; the effective angle of 
attack is given by equation (28). 
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